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We introduce a search method for a new class of gravitational-wave signals, namely long-duration 
O (hours — weeks) transients from spinning neutron stars. We discuss the astrophysical motivation 
from glitch relaxation models and we derive a rough estimate for the maximal expected signal 
strength based on the superfluid excess rotational energy. The transient signal model considered here 
extends the traditional class of infinite-duration continuous-wave signals by a finite start-time and 
duration. We derive a multi-detector Bayes factor for these signals in Gaussian noise using J^-statistic 
amplitude priors, which simplifies the detection statistic and allows for an efficient implementation. 
We consider both a fully coherent statistic, which is computationally limited to directed searches 
for known pulsars, and a cheaper semi-coherent variant, suitable for wide parameter-space searches 
for transients from unknown neutron stars. We have tested our method by Monte-Carlo simulation, 
and we find that it outperforms orthodox maximum-likelihood approaches both in sensitivity and 
in parameter-estimation quality. 



I. INTRODUCTION 

Gravitational waves (GWs), predicted by Einstein's 
General Theory of Relativity, are expected to be emitted 
by spinning neutron stars (NSs) with non-axisymmetric 
deformations or currents. Traditionally we distinguish 
two categories of GW emission from NSs: 

(i) short "burst-like" GWs from NS oscillations, e.g. 
fundamental (f-), pressure (p-), Rossby (r-) or w-modes, 
which could be excited by a NS glitch. These signals 
would typically be in the kHz frequency range and would 
be damped on timescales of milliseconds (see [HE])- Vari- 
ous data-analysis methods for such high-frequency bursts 
have been developed (e.g. see [2111]) and recently a search 
for a GW burst from a glitch in the Vela pulsar has been 
performed on LIGO data [5]. 

(ii) long-duration "continuous waves" (GWs) from non- 
axisymmetric deformations or currents in spinning NSs. 
These CWs are quasi-sinusoidal with a well-defined, 
slowly varying frequency /, which is typically of order 
of the NS spin frequency v. in particular f = 2v for 
"mountains" or precession, additionally f i/ for pre- 
cessing NSs, while / ~ 4i//3 for r-mode oscillations. In 
the past decade a number of data-analysis methods have 
been developed and applied to searches for GW signals in 
the data of ground-based detectors. These searches typi- 
cally come in two flavors, either fully coherently targeting 
known pulsars at f — 2iy, or semi-coherently searching for 
CWs from unknown NSs in a large parameter space of 
frequencies and sky-positions. See [6j for a review of the 
astrophysical models, data-analysis methods and results 
of searches for GWs, and for further references. 

The traditional GW model assumes that these sig- 
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nals are truly continuous in the sense of a quasi-infinite 
duration, or at least of longer duration than the avail- 
able observation time, typically Tobs ^ 1 ^ 2 yr. This 
postulates a spinning neutron star with a quasi-stable 
non-axisymmetric deformation, such that the relaxation 
to the axisymmetric thermodynamic equilibrium hap- 
pens on long timescales ^ Tobs- Alternatively, the non- 
axisymmetry can be driven by external influences such as 
accretion in low-mass X-ray binary systems (LMXBs), or 
by GW-driven instabilities such as unstable r-mode oscil- 
lations. Both are complex dynamical processes that need 
to be perfectly stable in order to produce a traditional 
CW signal. 

In the wide gap between burst-like O (ms) and truly 
continuous O (oo) signals, spinning NSs can reasonably 
be assumed to emit "GW" signals of intermediate dura- 
tion of O (hours — weeks) . We refer to this third category 
of GWs from spinning NSs by the oxymoron "transient 
CWs" . We can give three plausibility arguments for why 
it would be worthwhile to explore this new parameter 
space: (a) Young NSs exhibit enigmatic sudden spin-up 
events called "glitches", followed by a relaxation phase 
with timescale between days to months. This shows that 
internal dynamical processes on these timescales do ex- 
ist in NSs. (b) Equilibrium NS configurations are ax- 
isymmetric and GWs from quasi-stable deviations are 
expected to be weak: deformations strong enough to 
produce detectable GWs on Earth would therefore seem 
more likely to be associated with "catastrophic" tran- 
sient events, (c) Our astrophysical understanding of the 
universe is incomplete. Independently of any astrophys- 
ical models, this range of parameter space is currently 
not covered by any other searches and should therefore 
be investigated. 

The transient-CW signal model is very similar to the 
standard CW model, namely quasi-sinusoidal emission 
with a well-defined, slowly varying frequency of the order 
of the NS spin frequency, / ~ O (z^). Contrary to GWs, 
however, the signal has a definite start-time i° and a 
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finite efi^ective lifetime Tg, and the signal amplitude ho 
can be modulated by a window-function, for example an 
exponential decay. 

While the present study is concerned with extend- 
ing the traditional O (oo) CW search to finite durations, 
there is an independent and complementary effort under- 
way (called "STAMP" ) to extend the traditional O (ms) 
"burst" searches to longer durations [?]• This "long 
burst" method does not assume a parametrized signal 
model, but instead tries to find connected time-frequency 
patterns in the cross-correlation power between multiple 
detectors. These transients can span any duration from 
seconds to weeks. This approach will be computation- 
ally cheaper and more robust towards a wider class of 
unmodeled GW signals, but would therefore also be ex- 
pected to be less sensitive towards the particular class of 
transient- CW signals considered here. 

The plan of this paper is as follows: in Sec.[ll]we discuss 
the astrophysical motivation for transient CWs, includ- 
ing a simple signal-to-noise ratio (SNR) estimate based 
on the superfiuid excess rotation energy in NSs. Sec. |III| 
introduces the parametrized transient-CW signal model, 
and Sec. |IV| develops the coherent and semi-coherent 
Bayesian search methods for these signals. Sec. |V] 
presents numerical Monte-Carlo results on the detection 



efficiency, and Sec. VI illustrates the performance on pa- 
rameter estimation. Sec. |VII| gives a concluding sum- 
mary. Details on implementation and computing cost 
estimates are given in the appendix \K\ 



II. ASTROPHYSICAL MOTIVATION FOR 
TRANSIENT CWS 

Concrete astrophysical predications 

There are currently very few concrete predictions in 
the literature regarding the existence and properties of 
potential transient-CW signals from spinning NSs. No- 
table exceptions are the recent studies [SI IS] of GW emis- 
sion from non-axisymmetric Ekman fiow during the post- 
glitch relaxation phase, which typically lasts for days to 
months (e.g. see [IH])- The authors conclude that GWs 
from this mechanism could be detectable with second- 
or third-generation ground-based detectors. Another in- 
teresting recent idea [TT] suggests that giant magnetic 
flares in magnetars could trigger polar Alfven oscillations, 
emitting GWs at around ^ 100 Hz and lasting for days 
to months, although the details of the underlying physics 
are uncertain [12| . 

One can further speculate on a number of potential 
transient-CW emission mechanisms from spinning NSs. 
Many "classical" CW mechanisms discussed in the litera- 
ture (e.g. see [5] for references) typically have large uncer- 
tainties on the lifetime of the emission, and are therefore 
potential sources of transient-CW signals. For example, 
free precession of NSs occurs when the spin-axis is mis- 
aligned with the axis of symmetry, and has long been con- 



sidered a possible mechanism for CW emission (e.g. |13j). 
A more detailed analysis [M] concluded that the emission 
from free precession could be damped on a (highly un- 
certain) timescale of a few weeks to years. Young NSs 
could be deformed by extreme magnetic fields and re- 
sult in stronger transient CWs from damped free pre- 
cession |15| . Newly-born magnetars with strong toroidal 
fields could be subject to a magnetic instability produc- 
ing strong GWs on the timescale of several days [IB] . 
Similarly, the strength and timescale of the r-mode GW 
instability remain highly uncertain (e.g. see |17| for a 
review), despite a number of studies over more than a 
decade. This instability could operate on timescales of 
days to months in newly-born NS. As an example, see 
'18' for recent work on the large variety of possible sce- 
narios for the r-mode instability and spindown-evolution 
taking into account nonlinear mode couplings. 



Energetics and SNR of transient CWs 

It is instructive to consider the general relation be- 
tween GW energy emitted and the average expected SNR 
for transient CWs. This can be useful to derive an 
order-of-magnitude estimate of the expected transient- 
CW SNR in a simple toy-model. 

In the following estimate we assume that the CWs are 
emitted by a non-axisymmetric deformation e(t) of the 
quadrupole moment, emitting GWs at frequency /. See 
also [13 for a discussion of the relevant relations in the 
case of r-modes. The total energy emitted in GWs is 
^GW = / -^GW dt, in terms of the GW luminosity Lew 
(e.g. see [B]), which is 
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(1) 



where / is the axial moment of inertia, and e{t) is the 
dimensionless deviation from axisymmetry of the spin- 
ning NS. G is Newton's gravitational constant, and c is 
the speed of light. We can write the corresponding signal 
amplitude /io(i) at the observer as 



ho{t) = —i ^<t), 



(2) 



where d is the distance to the NS. Combining Eqs. ([T]) 
and ([2|, one can eliminate e{t) and write the total GW 
energy emitted during a time span T as 



27^2 c3 

Egw = ~^~Q~ f 



hl{t)dt, 



(3) 



assuming a roughly constant average frequency /. The 
expected optimal signal-to-noise ratio (SNR), which will 
be defined in Eq. ( 42 ) , can be averaged over sky-position 



and signal polarizations [13 [20] ; to yield 



{pD 



25 S{f) 



hl{t) dt, 



(4) 
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where S{f) is the (single-sided) noise power spectral 
density (at the signal frequency /) of the detector (or 
multi-detector combination, see [501 HI] ) • By combining 
Eqs. ([3| and Q, we obtain the average optimal SNR of 
a transient CW in terms of the total emitted GW energy 
i?GW, namely 



Egw 



2G 



57r2 c3 S{f) p dP 



(5) 



which agrees with the result given in [52]. An analogous 
expression was derived in 23 for the case of r-modes. 
This expression is interesting for two reasons: (i) at fixed 
frequency /, the expected SNR only depends on the to- 
tal GW energy emitted, and not on the timescale of the 
emission, i.e. short-strong or long-weak transient CWs of 
the same total energy typically result in the same SNR, 
and (ii) at fixed emitted -EgWj the optimal SNR po de- 
creases linearly with increasing GW frequency /, i.e. for 
the same transient GW energy, transients from slowly 
rotating NS would be easier to detect than from fast ro- 
tators. 

For convenience we Refine the dimensionless root- 
mean-square amplitude /iq over the timescale T as 



hl(t) dt. 



(6) 



Note that this quantity is not to be confused with the 
"root-sum-squared" amplitude h-^ss often used to char- 
acterize signal strength in the context of burst searches 
(e.g. see [3]), which is defined as /i^g^ = / h^{t)dt. The 
difference is that h^ss has dimension ^/t and refers to the 



measured strain h{t) in a given detector (see Eq. (16l), 
which is rapidly oscillating at frequency /. The dimen- 
sionless ft-o, on the other hand, refers to the intrinsic sig- 
nal amplitude ho{t), which for a transient CW would be 
slowly varying on a timescale T. From the definition ( 42 ) 
of the (optimal) SNR, one can see that = {2/S)hf^, 
and therefore Eq. Q also yields the relation 



(7) 



The quantity /i§ T is proportional to the GW energy ^ , 
as well as the average optimal SNR Q, namely 



hlT (X (pI) S cx 



GW 



(8) 



Transient CWs from superfluid excess energy 

As a simple toy-model for the energy available in prin- 
ciple for transient CWs from a spinning NS, we consider 
the standard two-fluid model, which is at the core of 
current attempts to understand pulsar glitches. See [10] 
for an overview of the phenomenology of observed pulsar 
glitches and the basic two-fluid model. 



The observed pulses (arriving with frequency v) are 
commonly associated with the rotating NS magnetic 
field, which is attached to the crust and the normal fluid 
interior, both rotating with angular velocity Q — 2ttv. 
This normal component has a moment of inertia /c, be- 
lieved to form the bulk of the total moment of inertia, 
I ~ Ic- The basic two-fluid model assumes that some of 
the interior neutrons are superfluid, forming an indepen- 
dent component that rotates at an (unobserved) angular 
velocity Us, and which has a moment of inertia Is, typi- 
cally believed to be of order Is/Ic ^ 10^^ (e.g. see [TU]). 
The normal NS components slow down at an observed 
rate O = 27ri' due to losses of angular momentum from 
the electromagnetic emission or interactions with the sur- 
rounding medium. The superfluid, on the other hand, is 
believed to be weakly coupled to the normal component 
(therefore ils ~ 0) and would continue to spin with an- 
gular velocity ils, until the "lag" Afl = fig — = — f2 At 
between the two components reaches a critical level (see 
[M] for a more detailed study including superfluid cou- 
pling effects) . The timescale At here could correspond to 
the inter-glitch period, if one assumes that every glitch 
restores perfect co-rotation between the fluids. At this 
"breaking point" some type of instability is believed to 
occur, resulting in the transfer of angular momentum 
from the superfluid to the normal fluid. This would pro- 
duce the visible "glitch" , i.e. a sudden spin-up SQ — 2t:5v 
of the observed pulse frequency v. The details of this in- 
stability are poorly understood and various models have 
been suggested in the literature, such as the crust break- 
ing due to the strain exerted by pinned vortices (e.g. see 
[251 )■ the vortex array becoming unpinned due to the 
Magnus force (e.g. see [26]), or a two-stream instability 
developing due to the dynamical coupling of the two flu- 
ids [571 HH] • See also [53] for an alternative mechanism to 
vortex pinning that can produce crust strain. Similarly 
unclear are the physical mechanisms involved in the ensu- 
ing "relaxation" back to a state of pre-glitch steady-state 
spin-down, which typically occurs on timescales of days 
to months (see p^). 

During a glitch, the observed normal fluid would spin 
up by 5il, while the superfluid would spin down by 5ils, 
such that angular momentum is conserved. Following [24] 
and assuming that the moments of inertia are constant 
during a glitch, we have 



(9) 



where typically the observed spin-up is up to 5il/il ^ 
10"^. Note that for flducial values of Is/Ic 10"^ 
this implies that the superfluid spins down by 5fls/il ^ 
— 10^* during a glitch. Besides this angular-momentum 
transfer in the glitch, there is some excess energy iS'giitch 
left. Assuming that co-rotation between the two fluids is 
restored after the glitch, namely 6ils = —Ail -\- 5il, the 
glitch excess energy can be shown to be 



E, 



glitch 



^Is6nl-^^Ic6il'^^Is6ill, (10) 
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which agrees with the result in [21], where in the last 
step we assumed the fiducial scales Is/Ic ~ 10~^ and 
therefore 5^/5Q.s ~ 10~^, which makes the second term 
negligible. Using Eq. ([o]), this can also be expressed as 
£^giitch ~ ^/c^rJAfi. This energy would have to be dis- 
sipated directly in the glitch, for example by exciting 
oscillation modes, which would radiate GWs on short 
timescales, as considered in [31 [5], or by heating up the 
NS (e.g. see [3D]). 

However, after a glitch the NS generally relaxes back 
to a steady-state spindown on a timescale of days to 
months, and therefore some (if not all) of the initial ex- 
cess energy stored in the faster-rotating superfluid could 
be driving GW emission on this relaxation timescale of 
O (days — months). The recent studies of post-glitch 
spin-up of the fluid NS core by non-axisymmetric Ek- 
man flow [HI |9] provide one concrete example for exactly 
such a mechanism. It is also conceivable that this hidden 
energy in the superfluid "flywheel" can lead to transient 
CW emission directly, via an internal instability, without 
transferring its angular momentum to the curst first, i.e. 
without a visible pulsar glitch. 

The well-known spindown upper limit on CW emission 
from known pulsars ^ assumes that the total steady- 
state spindown energy of the pulsar is emitted as GWs. 
Similarly we can consider a "superfluid upper limit" on 
transient CWs, where the total superfluid excess energy is 
converted into transient CWs, for example by sustaining 
some non-axisymmetric process that leads to emission at 
frequencies / ~ 0{v). The excess rotational energy 
of the superfluid is 



E, = -h{ni-n')^An^hvi\v, 



(11) 



where in the last expression we dropped the second-order 
energy term (10), which is smaller by Ail/f2 w 10"''. 



In the following we return to using the spin frequency 
V instead of the angular velocity VI = 2-11 v. Note that 
this expression happens to be numerically very simi- 
lar to the "starquake" glitch-energy derived under the 
assumption of a spin-up caused by a reduction in the 
crust's moment of inertia (e.g. see [5]), which results in 
-E-quake = liT^ I V 5v , and substituting 5v — (/s//)Ai/ 
from (|9]), we obtain 2i?quako = Es. For fiducial values 
of / ^ lO^^kgm^, Vela spin frequency v ^ 10 Hz and a 
large ghtch oi 5v/i> 10~^, this yields Eg ^ Ax 10''^ erg. 

We have no direct observational evidence about the 
size of the lag Ail = 2ti/S.v or the superfluid angular 
velocity fig. In the simplest models one typically assumes 
the lag to be reset to zero after every glitch, starting to 
build up again due to spindown v. However, this is not 
necessarily the case and the lag could accumulate over 
longer timescales and be correspondingly larger. 

To obtain an upper limit estimate on the signal 
strength, we can equate E'gw of Eq. (|3| with E^ given in 
Eq. (Ill, assuming that the GW is emitted at a frequency 
f = 2i> (corresponding to a non-axisymmetric deforma- 
tion). This corresponds to the extreme case where the 



built-up excess superfluid rotational energy drives emis- 
sion of a transient CW. Combining this with Eq. ([6]), we 
can obtain a superfluid transient-CW upper limit esti- 
mate in the form 



1 



(12) 



Alternatively, if we assume that the total excess angular 
momentum is transferred in a glitch, we can use Eq. ^ 
to rewrite this in terms of more directly observed glitch 
quantities, namely 



(13) 



where / is the axial NS moment of inertia, and 6v/v is 
the observed glitch spin-up. As a third alternative, we 
can use the simple two-fluid model for the buildup of the 
lag, namely Av = {vlAt, and parametrize the unknown 
superfluid moment of inertia via Ig = "f I (with fiducial 
value 7 ^ 10^^), and obtain the relation 



(14) 



in terms of the well-known "spindown limit" /igj, given 
by 



(15) 



which is derived from the assumption that the NS spin- 
down energy is completely converted into GWs. The re- 
lation (14) shows that the superfluid transicnt-CW up- 



per limit can be similar or even higher than the usual 
pulsar spindown limit: a fraction 7 of the steady-state 
spindown energy is accumulated in the superfluid over 
an inter-glitch timescale At, and is released on a short 
timescale T. Assuming 7 ~ 10~^, inter-glitch periods of 
At ^ 1 yr and transient-CW timescale of T ~ 2 weeks, we 
see from Eq. (14 1 that for these values Hq « hsd- This is 
interesting as it indicates that the young (and glitching) 
pulsars with the most promising spindown upper limits 
on CW emission might also be the most attractive tar- 
gets for directed transient-CW searches, for example the 
Crab pulsar. Vela, and J0537-69 (e.g. see Fig. 4 in [6]). 



III. SIGNAL MODEL FOR TRANSIENT CWS 

The family of transient-CW signals considered here is a 
straightforward generalization of the traditional infinite- 
duration CW model ([13, 31 , allowing for a finite du- 
ration and non-trivial time-evolution of the overall am- 
plitude ho. The set T of extra transient-CW parame- 
ters therefore consists of the start-time t", characteris- 
tic timescale r, and the type of window function w, i.e. 
T = {tu, r}. The transient-CW signal family therefore 



5 



simply consists of a window-function gT^{t;t'^ ,t) applied 
to the standard CW signal model, namely 



/i^ it; A, A, T) - g^it; r) A'' [t; A) 



(16) 



were we use implicit summation over the four amplitudes, 
fi — 1, ... 4, and X is an index over different detectors. 
The four canonical amplitudes A'^ are functions of the 
CW amplitude ho, polarization angles cosi, V' and the 
initial phase ^o, i.e. A^ = A'^{ho,cosL,tp,4>o). The cor- 
responding basis functions {t; X) are found, for exam- 
ple, in |32j . but their explicit form is not relevant to our 
discussion here. The set of Doppler parameters A deter- 
mines the time evolution of the signal phase, for example 
the source sky-position n and the GW frequency f{t), 
which is generally allowed to be slowly varying with time. 
If the CW source is a neutron star in a binary system, A 
would also include the relevant orbital parameters of the 
system. 

In the following we restrict our attention to two simple 
types of transient window functions (t) , namely either 
rectangular, denoted as = r, i.e. 



?r(i;t°,r)^| 



_ / 1 if t e [t",to + r] 

otherwise , 



(17) 



or exponentially decaying, denoted as tu = e, namely 

g.it;t",r)^lf'-''''' ifte[^",^o + 3.] (^3) 
' ' ^ \ otherwise, ^ ' 

where we somewhat arbitrarily truncated the exponen- 
tial window at an e-folding of 3, in order to simplify the 
practical implementation of this window. This trunca- 
tion gives the window a finite duration of 3r, and at the 
truncation-point the amplitude has decreased by more 
than 95% and we can neglect the corresponding loss of 
SNR. 

Following the notation of [21] |33j , we use boldface to 
denote multi-detector vectors, i.e. {x}-'^ = denotes 
the set of data-streams x(t) from different detectors X. 
We can now conveniently absorb the window-function 



g(t}) in Eq. ( 16 1 into the definition of new transient basis 
functions h! , namely 



h!At;\,r)~g^{t;t\T) h^{t;\) . 



(19) 



If we denote 9 the set of all signal parameters of our 
search, i.e. 



= {A,\,T}, 



(20) 



then we can write the transient signal model ( 16 1 now 
more compactly as 



h{t;d)^A^h!^{t;\,r). 



(21) 



IV. DETECTION METHOD: ODDS RATIO 

A. Hypothesis testing framework 

Based on observed data x, we want to decide between 
two hypotheses: under the noise hypothesis the ob- 
served data consists only of Gaussian stationary noise n, 
and under the signal hypothesis the data contains 
in addition a transient-CW signal h{t;d) of Eq. (21 1, 
namely 



-Hg : x{t) = n{t) , 
Hs ■■ x{t) = n{t) 



h(t; 9) for any 9 e 



(22) 
(23) 



where P denotes the signal parameter space. Note that 
the signal hypothesis ( 23 1 is incomplete without the spec- 



ification of a probability distribution for the unknown 
signal parameters 9 over their parameter space P, i.e. a 
prior probability ^(^17^3). Note that for simplicity our 
notation does not explicitly distinguish between proper 
probabilities and probability densities. This difference 
is implicit in the type of argument, namely whether it is 
discrete, such that J2i P [xi\T) — 1, or continuous, where 
J P {x\X) dx = 1. Furthermore, we sometimes (but not 
always) explicitly state I as a conditional in probability 
statements, denoting the full set of remaining relevant 
prior model assumptions that the probability statement 
depends on. 



B. Gaussian noise and scalar product 



For quasi-sinusoidal CWs (161 in stationary Gaussian 



noise we can define a multi-detector scalar product 
using the notation of [33], as 



{x\y)^2Y,S-^\f) 

X 



x''{t)y''{t)dt, (24) 



where Sx{f) is the (stationary) single-sided noise power 
spectral density in detector X, which is assumed con- 
stant over a narrow frequency band around the signal 
frequency /. This allows us to write the likelihood for 
the data x in the Gaussian noise-case ( 22 ) as 



P{x\Hc 



-i(x|x) 



(25) 



where k is a normalization constant. In the presence of a 
signal h{t; 9) with parameters 9, subtracting this signal 
from the data x results again in pure Gaussian noise n, 
i.e. n = X — h{9), and therefore 



Pix\ns,( 



(26) 



The likelihood for the data x containing any signal h{t; 9) 
with 9 G V drawn from the prior P (6'|'Hs) can easily be 
obtained (e.g. see (Ml) as 



Pix\Hs)^ I P{x\Hs,9) Pi9\Hs) d9 , 
Jp 



(27) 
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which is often referred to as the marginal likelihood (and 
sometimes as evidence). 



C. Odds ratio and Bayes factor 

We can express the odds ratio OsG between signal- and 
Gaussian-noise hypothesis for the observed data x as 

OsG\x)= I X ^d<s.g(x) (28) 



p{nG\x) 



p{nG\i) ' 



where we used Bayes' theorem, namely P (a|fo) P (6|I) = 
P (6|a) P (a|I), to express 0'qq,(x) as a product of the 
prior hypothesis odds and the so-called Bayes factor 
Bsg{x), defined as 



Bsg{x) 



P{x\'Hg) 



f c{x;e)p{e\ns) d0, (29) 

Jp 



where we used Eqs. (25), (26) and defined the standard 
likelihood ratio as 



_ Pix\ns,6 

p{x\nG) 



C{x; 9) = ^ 7'- -'"/ ^ ^(.\H0))-UH0)Me)) ^g^) 



Inserting the transient-CW signal model (21), we can 
write this as 

\ogC{x-e)=A^^x'^~\A^M'^,A\ (31) 



where we defined 



x'^{\T)^{x\h'^) 



(32) 
(33) 



generalizing the corresponding CW quantities, e.g. see 
|34j . to the transient signal model. 



D. Maximum-likelihood: the J-"-statistic 



Contrary to the marginalization in Eq. ( 29 1 over un- 



known parameters 0, which follows from the axioms 
of probability, the orthodox "maximum-likelihood" ap- 
proach consists of an ad-hoc maximization of the like- 
lihood ratio C over the unknown parameters 9, i.e. we 
define the maximum-likelihood statistic as 



'CmlI^;) = max£(a;;( 



(34) 



transient-CW case, the transient J^-statistic is obtained 
explicitly as 



2F{x-\T) = x'M 



(36) 



where ^A'^^'^ is defined as the inverse matrix of A^^^ of 



Eq. (33) 



If the data x(t) contains a signal s{t]9s), such that 
X = n + s{9s), then we can write Eq. (32 1 as 



(37) 



and 



(38) 



with the obvious definitions n'^{X,T) = (n|/ij^) 

s'^{9,;X,T)^{si9,)\h'^{X,T)) , 

which depends both on the signal parameters = 
{^s)As,7^} and the matched-filter parameters {A,T} of 
the "template". 

Gaussian detector noise n{t) has zero mean, i.e. 
E [n] = 0, and therefore E [n^] = and E [x^J = s'^. 
One can also show that the corresponding covariance is 
E [n'^nl] = M'^^{X,T), and therefore 



E[x'^xl] =M'^, 



(39) 



Using this together with Eq. (36) one can further show 
that 2T follows a -distribution with four degrees of 
freedom and non-centrality parameter p^, i.e. 

E [2 J-] = 4 + (40) 

where the signal-to-noisc ratio (SNR) p is expressible as 



(41) 



We see from Eq. (|38|) that the SNR will depend in a 
complicated way on the offset between signal parame- 
ters 6*8 and template parameters {A,T} (see [33J for the 
non-transient CW case). In the special case of perfectly 
matched template parameters, i.e. A = As and T = 
we obtain the so-called "optimal SNR" po, which can be 
expressed as 



p2(0,) = ^A^x;^^(A„rs) Ai = {8\s) 



E. Choice of signal priors 



(42) 



In order to fully define the Bayes factor ( 29 ) , we need 



to provide a complete signal hypothesis including the 
prior probabilities P (^jHs) for the signal parameters 

= {AA,r}. 



Given the explicitly quadratic dependency on A^ in 



Eq. (31 1, this maximization can be performed explicitly, 
which 



1 results in 

In Lmi.{x) 



max J-(x \ X.T) . 
{A,r} 



(35) 



were we encounter the well-known "J^-statistic" , which 
was first derived in fT3] for CW signals. In the present 



1. Prior on Doppler-parameters A 

For simplicity we assume that the Doppler-parameters 
A are independent of amplitude- and transient parame- 
ters {A, T}, so we can factor the full parameter prior 
P{9\Hs) into 



P((?|Hs) = P(A|Hs) P{A,T\'Hs) , 



(43) 



7 



and so the Bayes factor ( 29 ) now reads as 

Bsg{x) = J Bsoix; A) P (A|Hs) d"A , (44) 

in terms of a "targeted" Bayes factor Bsg{x; A) for a 
single Doppler point A, namely 

Bsg{x; A)= J C{x:A,T,\)PiA,T\ns) dAdT. (45) 

In the following we will focus exclusively on the targeted 
Bayes factor, and sometimes drop A for simplicity of no- 
tation. The generalization to parameter searches over A 



is straightforward as given by Eq. ( 44 ) 



2. Prior on transient parameters T 

Astrophysically it would make sense to assume that the 
amplitude hg of a transient CW is related to its timescale 
T. For example, it might be reasonable to suspect that 
stronger transient GWs have shorter duration and vice- 
versa, according to some prior on the total transient-CW 
energy emitted and on the distance of such sources, cf. 
Sec.jn] However, given the current lack of concrete astro- 
physical predictions to base such priors on, we assume a 
simple independent prior P (Tj'Hs) for the transient sig- 
nal parameters T. A naturally simple choice consists 
of independent uniform priors within some appropriate 
time- windows, i.e. 



P(t"|Hs) =[/(Ci„,Cin + At°) 

P{T\ns) = U (r^in, T,ni„ + Ar) 



(46) 



where we write uniform probability densities as 
U {a, b) = l/{b — a) for the parameter falling inside [a, b] 
and zero otherwise. We would also need priors on the 
window type, e.g. P — r|Hs) = -P (ro = ej^s) = in 
order to marginalize over w. However, for simplicity we 
will often assume a particular window-type vj as given 
and only marginalize over t'^,T. The effect of assuming 
an incorrect window-function within tu e {r, e} is numer- 



ically studied in Sec. V B and appears to entail only mild 
losses of detection power. 



3. Prior on amplitude-parameters A 

Physically reasonable priors on the angle parameters 
{cos 6, 0o} are relatively easy to obtain if we assume 
ignorance about the orientation of the spinning neutron 
star (e.g. see [33] for a more detailed discussion), namely 
by symmetry one can obtain 



P{cosL\ns) = U{-l, 1) , 

Pirns) {-II) 

P(0o|Hs) = f/(0, 2tt) . 



(47) 



The choice of prior for the overall amplitude parame- 
ter ho is less obvious: one could choose a scale-free Jef- 
frey's prior, or a simpler uniform prior on some physically 
meaningful domain. 



The downside of the isotropic amplitude priors (47) 



is that the resulting marginalization over A cannot be 
performed analytically, and computing the Bayes factor 
would require a numerical integration over A in every 
point A,T, which will be computationally prohibitive. 
However, as shown in [34j . by using an unphysical uni- 
form prior on the A-vector A'^, one can analytically 
marginalize over A and obtain a Bayes factor ( 45 1 ex- 
pressed in terms of the well-known J^-statistic ( 36 ) . 



This choice has the major advantage of simplicity and 
computational efhciency, while incurring only small losses 
in detection power compared to the more physical prior 
(47), as shown in [34]. Given that fast and mature codes 
exist to compute T on real detector data (cf. f20l'35]), we 
choose this prior as a convenient practical approximation. 

However, as will be seen in the following, the original 
formulation of the uniform ^'^-prior in [34] leads to a 
somewhat weak detection statistic for transient-CW sig- 
nals and can be improved by a minor modification. The 
original "constant-Zimax" prior was defined as 



PiAf^lHi 



C 




if tio{A) < K 
otherwise , 



(48) 



where /imax is a maximum cutoff amplitude needed in 
order to normalize the prior. 



The ^-integration in Eq. ( 45 ) can be performed analyt 



ically with this prior, provided the data does not cause 
the likelihood to peak close to the upper cutoff /imax- 
Namely, if the value of the integrand C{x]A) is already 
negligible at the cutoff boundary, we can extend the do- 
main to infinity and obtain a 4-dimensional Gaussian in- 
tegral, namely 



C{x;A,T)Cd'^A: 



^/\M\ 



(49) 



where = detA^ is the determinant of the matrix 



Mfj,i, of Eq. (33). Using this approximation, we can 



therefore write the Bayes factor ( 45 ) as 



(50) 



The Jacobian J of the coordinate transformation = 
A'^{hQ, cost, -0, 0o) is J = j/iq ~ (s^^ Elj), and 

therefore dt4 = J dho dcos t dip d(j)o ■ We can now deter- 
mine the normalization constant C as 



1 = J P{A\ns) d^A 



C- 



27r2 
35 



(51) 



and obtain the constant-Zimax Bayes factor explicitly as 
70 



At" At 
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where we assumed a fixed window type w. The antenna- 
pattern weighting factor |A^|~-'^/^ generally depends on 
the sky-position n, the transient-window type tu, start- 
time i° and timescale r, as seen from Eq. (33). In the 



case of a fully targeted search with fixed v\t and A, 
as discussed in [34;, the weighting factor is constant and 
plays no role for the power of the detection statistic. The 
fully targeted Bayes factor is therefore strictly equivalent 
to the detection power of the J^-statistic. 

Interestingly, in the case of transient signals we find 
that the presence of this antenna-pattern weighting fac- 
tor in Eq. ( 52 ) seems to degrade, the detection power, 
and for some choices of parameter ranges i^scl/i per- 
forms worse than the maximum- likelihood statistic (35), 
namely 

■^max{x; vj) = max J-{x] m,t'^,T) . (53) 

{tO,r} 

However, a simple modification of the cutoff boundary 
in Eq. ( 48 ) allows us to eliminate the antenna-pattern 
weighting factor in Eq. ( [50| . Namely, by intro- 

ducing an "SNR-scale" p = ho and using a cutoff 

P < Pmax as the outer boundary of the domain instead 
of the amplitude-cutoff /iq < /imax, the modified ad-hoc 



"J^-statistic prior" is now 



P{A^\Hs,T) 



C 




if p{A) < p„ 
otherwise , 



(54) 



which results in the new transient Bayes factor 
70 



BAx; A, tn) = ■ / e^(-^^^^) dt° dr , (55) 

Pmax '-^^ 



-Bsclp^^^^as our 



where in the following we denote Bjr 
detection statistic of choice. 

As discussed in Sec. |V A[ numerical simulations show 
that this detection statistic is more powerful than both 



as well as the orthodox maximum-likelihood 



statistic. More work would be required to study these 
priors in more detail, and to develop a more physical 
choice of priors that would be equally practical. 



where we have omitted multi-detector summation for 
simplicity of notation. It will be difficult to make fur- 
ther analytic progress with this expression in the general 
case, but we can analyze the interesting special case of 
the rectangular window function (17), i.e. = -nj = r. 



In this case the window-functions simply truncate the 
integral, and so we obtain 



2^ 
S 



Kit) Kit)dt, 



(58) 



to 



where we defined to = max(tg,t''), and ti = 
min ((tg -I- Ts), (i° -f r)). Note that [^0,^1] denotes the 
rectangular overlap between % and T, and in the above 
integral we assumed ti > to, otherwise the expression is 
zero. 

In order to simplify this even further, we note that the 
antenna-pattern matrix in Eq. ( 33 ) can be written as 



M'. 



2r 



(59) 



where we defined m^^^ ~ (/i^ /ii/)s, in terms of the multi- 
detector time-average introduced in [Eq. (59)] in 
|33) . The antenna-pattern functions are periodic with 
period of a sidereal day, and so the average rn^i/ will 
be weakly oscillatory and converges to a constant for 
T ^ Id. Let us therefore approximate m^^ as constant 
for fixed A, i.e. m^i, « w^i,, which allows us to write 



Eq. (58) as 



with TA = [ii — <o] 



(60) 



where [...]'*' is the positivity operator, defined as [a;]^ = 
X for a; > and zero otherwise. Therefore ta is the 
length of overlap between signal and template windows. 
Using M!^^y « {2t / S)m^v, we obtain the mismatched 
SNR after substituting into Eq. (41), namely 



^ 2 



(61) 



F. SNR-loss due to rectangular-window offsets 

Let us consider the effect of an offset in transient- 
window parameters 7" from the signal parameters 7^, as- 
suming perfectly-matched Doppler parameters, i.e. A — 
Ag. The dependence of the matched-filter SNR (41) on 
Doppler offsets As — A signals has already been studied 
in great detail, e.g. see [STl [551 IBS] . 



In the following we drop A and write Eq. ( 38 ) more 
explicitly as 



s'^{A,%,T) = a: {K{%)\h'^iT)) 



(56) 



and using Eqs. ( 19 ) and ( |24[ ), we can further expand this 

as 

" ' K{t)h^{t)g^Xt;t'i,T,)g^{t-t\T)dt, (57) 



S 



in terms of the perfectly-matched "optimal SNR" po, de- 
fined in Eq. (42). Note that always ta < min(T, Tg), and 
equality only holds in the perfect-match case. Expressing 
this in terms of the usual definition of mismatch m, we 
obtain 



;i°,r)^^ 



Pi 



1 - 



'A 



(62) 



The behavior of this approximate mismatch function 
and the corresponding measured SNR loss is illustrated 
in Fig. [1] for a start-time t^ = 5 days and duration 
T — h days |46| . Contrary to the well-known mismatch 
behavior in Doppler parameters A, the transient mis- 
match metric is not differentiable at the perfect-match 
point T = Tg = TA, as seen in Eq. (62 ) and Fig. [I] where 
the mismatch has a kink. Therefore we cannot Taylor- 
expand the mismatch around the signal location and 
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FIG. 1: Approximate mismatch (62 I (dashed lines) and mea- 
sured SNR loss (solid lines) as functions of the rectangular- 
window offsets {t°, r}. The top left plot shows mismatch con- 
tour lines at m = 0.90, 0.75, 0.60 respectively, and the side- 
panels show cross-sections of / pi at fixed f° — (right 
panel), and at fixed r = Ts (bottom panel). 



define a metric tensor from the second-order quadratic 
form. 



Furthermore, the mismatch ( 62 ) depends not only 
on the parameter offsets, but also on the actual value of 
the signal duration Tg. Therefore the mismatch function 
is not constant over the parameter space. However, it is 
interesting to note that the mismatch behaves close to 
linearly around the point of perfect match, and we can 



obtain the first-order variation of the mismatch ( 62 ) as 



(63) 



which shows that the mismatch increases twice as fast 
for offsets in start-time f° than for offsets in timescale t, 
which is also seen in Fig.[l] Also, as seen m Fig. [TJ the 
parameters and r are correlated, and the iso-mismatch 
curves close to the signal are straight lines with steepness 
\dT/dt^\ = —2, as seen from Eq. (63|. 

These properties are important for covering the 
transient-parameter space with a template bank, where 
one tries to use the smallest number of templates while 
guaranteeing a certain minimum match for any sig- 
nal within the template-space. For our present pur- 
pose it will be sufficient to make sure that the finite 
step-sizes in t^,T are fine enough to ensure a reason- 
able approximation to the integral ( 55 ) . The worst-case 



mismatch (62) occurs for the shortest timescale Tg, so 
for Tmin = 0.5 days, say, and a time-sampling in steps 
of TsFT = 1800 sec, the worst case mismatch will be 



G. Semi-coherent Bayes factor 

Increasing the coherent integration time (or in our 
case, the maximal timescale r) in wide parameter-space 
CW searches over unknown Doppler parameters A typ- 
ically results in a dramatic increase in computing cost. 
The reason is that the likelihood function F becomes in- 
creasingly finely structured over Doppler parameters A, 
such that more and more templates need to be sampled 
in order to cover the parameter space (see (HU |33j l36]). 
This feature severely limits the feasible coherence time 
to ~ O (days) for fully coherent searches over unknown 
Doppler parameters. 

The usual approach to this problem is to abandon 



fully coherent integration of Eq. ( 24 1 over the whole ob- 
servation time and instead adopt a "semi-coherent" ap- 
proach: this is typically achieved by relaxing the con- 
straint of a consistent signal phase over the whole life- 
time. Namely, one splits the observation time Tobs 
into N segments of duration AT, such that Tobs = 
N AT, and requires phase-coherence only over each seg- 
ment AT. The initial phase 0o is part of the set of 
four amplitude-parameters A, but for simplicity one re- 
laxes the consistency-constraints for all four amplitude- 
parameters across different data-segments a;(i) with 
i = 1 . . . N. This corresponds to replacing the four un- 
known amplitude parameters A'^ by x 4 unknown am- 
plitude parameters {A'^-s^}^i. Using the product rule 
for joint probabilities of independent events, namely 

P {{xi^i-j}\ ■ ■ ■) — YliLi P (^(i)l ■ ■ ■)' ''^'^ express the 
corresponding semi-coherent transient Bayes factor (29) 
as 

BlU^) ^ Jd^J dTP{X,T\ns) X 

N 

Y[J C{x(,y, ^(,) , A, T) P (^(.) \ns) . (64) 



Using the J^-statistic prior ( 54 ) , we obtain 

N 



1 



AiO At \p 



70 

max 



(65) 



where we defined the semi-coherent sum Sjv as 



N 



5jv(a;; A,r) = ^ J"(a;(,); A,r) . 



(66) 



The semi-coherent sum can be shown to require substan- 
tially fewer templates in parameter space for the same 
amount of data analyzed e.g. see [371 [38] . and can in many 
cases achieve a higher sensitivity for wide parameter- 
space CW searches at fixed computing cost. In or- 
der to maximize the available computing power, such 
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searches are performed on large computer clusters and 
on Einstein@Home [47] . a large public distributed com- 
puting project (e.g. see [52] for an Einstein@Home search 
on LIGO data). 



Therefore the semi-coherent Bayes factor ( 65 1 allows 



for an efficient and computationally feasible transient- 
CW search over a wide Doppler parameter space A of 
sources with unknown sky-position, frequency and spin- 
down. More work is required to fully develop and study 
this approach, as our present analysis is mostly focused 
on the coherent method. A semi-coherent search for tran- 
sient CWs would be very suitable to be performed on the 
EinstcinQHome computing platform. 



V. DETECTION EFFICIENCY 

In order to quantify the detection efficiency of differ- 
ent statistics we perform Monte-Carlo studies with sim- 
ulated Gaussian noise and injected signals with parame- 
ters drawn according to their priors. Note that in order 
to implement this in the most efficient way, we do not 
perform end-to-end simulations from generated artificial 
data processed through the whole pipeline. Rather, we 
more directly synthesize draws of the different statistics 
from their known distributions, or compute them from 
draws of intermediate data-products. This is possible 
because we know the distribution of intermediate data- 
products in the case of Gaussian noise. The details of 
this procedure are described in Appendix |A 4[ 

A useful method to compare different detection statis- 
tics is to compare their "receiver-operator-characteristic" 
(ROC), namely the expected detection probability pdct 
versus the false-alarm probability for some signal 
population. This follows the spirit of the Neyman- 
Pearson criterion, which compares the detection prob- 
ability of different statistics at fixed false-alarm proba- 
bility. 

Note, however, that the traditional way of plotting 
ROC curves, namely p^ot versus is somewhat unfor- 
tunate. We know a priori that any statistic of positive 
detection power satisfies pdot > Pfa, and pdot — Pfa cor- 
responds to a complete random guess. Therefore only 
the upper triangle of the plot is of any interest, and it 
is more informative to plot pdct — Pfa versus Pfa instead, 
which quantifies how much better any statistic performs 
compared to a random guess. 



A. Comparing different amplitude priors 

We first consider the detection efficiency of the two 
Bayes factors Bgclp^ 



of Eq. (55), and ^Sg|;j,„^^ of 
Eq. (|52|). These only differ by their cutoff boundary 



noise- and in the signal-case. We estimated the errors on 
Pdct using a jackknife estimator on 100 subsets, the re- 
sulting estimated la errors in the following ROC curves 
are always less than (i(pdot) ^ 0.02. 

The amplitude parameters are drawn according to 
their physical priors (47). For the first simulation we 
fixed the optimal SNR (42) to = 3. This is achieved 



by re-scaling Hq according to the resulting SNR for drawn 
values of {cosi, V', 0o}- We always use the same ranges 
for the search and for drawing signal parameters from, 
and we consider two sets of transient window parameter 
ranges, namely 

range I: r € [0.1,0.6] days, t° e [0,9] days, 

(67) 

range II: t e [0.5, 2.5] days, t° e [0, 6] days . 

The results of this Monte-Carlo simulation are shown 
in Fig. j2j and we see that the statistic Bjr = ^sclp 
seems to generally perform better than BsgI/j j as dis- 
cussed in Sec. IIVE3I For some choices of transient win- 
dow ranges (such as range I) , the latter actually performs 
worse than the orthodox maximum-likelihood statistic 
-^max, as seen in the upper plot in Fig. [2j We also note 
that the detection probability for signals of equal optimal 
SNR po is lower for range I than for range II. This can 
be understood from the substantially larger parameter 
space associated with range I, namely Tmax = 0.6 days in 
Tobs = 9 days, as compared to Tmax — 2.5 days in Tobs = 
6 days for range II. Therefore we can fit Tobs/T = 15 inde- 
pendent rectangular windows into the observation time 
in range /, while for range II this factor is only 2.4. This 
entails more independent trials and therefore a higher 
false- alarm probability for range I , which reduces the de- 
tection power for signals of the same SNR. 

In order to verify that these qualitative conclusions are 
not restricted to injections at constant SNR po, we also 
performed these Monte-Carlo simulations for injected sig- 
nals at constant amplitude ho / VS. Note that this results 
in a wide range of injected signal SNRs po due to the 
varying signal durations r. The results of these simula- 
tions are shown in Fig. [3] which qualitatively agree with 
Fig. [2] 

As expected, these results also show that the transient 
statistics are substantially more sensitive to transient- 
CW signals than a standard "infinite-duration" CW J-- 
statistic search over the full span Tobs- The quantitative 
advantage in recovered SNR depends on the details of the 
transient parameter space, as seen from the mismatch 
(62), namely setting t — Tobs and ta = Tg, we find the 



on the uniform-y^'* amplitude priors, as discussed in 
Sec. |IV E 3| For simplicity we consistently used a rectan- 
gular transient window (17) for both the injections and 
the search. We performed N = 10^ random draws in the 



recovered fraction of p^ is roughly proportional to the 
"duty cycle" Ts/T^hs of the transient signal with respect 
to the observation time Tobs- 



B. Comparing rectangular and exponential 
windows 



Another question of interest is how robust the detec- 
tion statistic B'^ is, which assumes a particular transient- 



11 




< 

I 

0) 
Pin 






0.6 




0.5 


< 






0.4 


1 


0.3 


0} 
T3 




0.2 




0.1 








range II 

ho = omvs 




Pf 



Pf 



FIG. 2: Detection efficiency of BsG|p„axi ^scUmax f^nd J-'max 
for injected signals at fixed optimal SNR of Po = 3, using a 
rectangular transient window (ots = r). The curve labelled 
'.Ftotai' refers to a standard CW J^-statistic search over the 
full data span. The upper plot corresponds to transient- 
window parameters drawn from range I , while the lower plot 
is for range II, see Eq. ( |67[ ). We only plot la error- bars for 
Bsclpniaxi which are representative for the size of the errors 
on the other curves. 



window type w, if the transient signal actually has a dif- 
ferent window type mg. We cannot answer this question 
in general, but it is instructive to study the simple case 
of type-mismatch between the rectangular and exponen- 
tial transient- window types. We inject signals with either 



rectangular (Eq. (17)) or exponential (Eq. (18l) transient 
window, and we perform the search using both rectangu- 
lar and exponential transient windows, respectively. The 
Monte-Carlo parameters are the same as in the previous 
section, and we only show the results for the transient pa- 
rameters range II of Eq. (|67| , the qualitative conclusions 
are the same for range I. The results of these Monte-Carlo 
simulations are shown in Fig. |4] We see that using the 
wrong window- type, i.e. w Wg, results in a loss of de- 
tection power in Bjr, as expected. However, the loss is 
quite moderate, which in practice would favor a search 
assuming the simpler and much more computationally 
efficient rectangular window type, zu — t (cf. Sec. A 3 1, 
if the search is computationally limited. If computing 
cost is not an issue, for example in a targeted search for 
known pulsars, one could perform both searches and then 



FIG. 3; Same as Fig. |2] with signals injected at fixed ampli- 
tude ho = 0.06\/5 and ho = 0.02\/S, respectively. 



marginalize over ro. More importantly, however, these re- 
sults suggest that we can hope to be reasonably sensitive 
also to transient signals with different time evolutions 
not contained in tu G {r,e}. 

It is interesting to note in Fig. |4] that in the case of 
injecting rectangular-window signals, i.e. Wg = r, the 
maximum-likelihood T^^.^ assuming an exponential win- 
dow, can outperform the more correct J^jJ^j^^-statistic. A 
similar effect is seen in the corresponding simulations for 
range I, which are not shown here. The origin of this 
"anomaly" can be traced to the fact that for the same 
timescale parameter r, the exponential-window wave- 



form ( 18 ) lasts three times longer than the rectangular- 



window waveform (17 1. Therefore, the parameter spaces 
range I and range II contain substantially more inde- 
pendent trials for the rectangular- window waveform com- 
pared to the exponential one. This results in lower false- 
alarm probabilities at fixed threshold for -F^^^^ compared 
to -Fjjjax- Although there is a loss in SNR due to window- 
type mismatch, this effect is partly weaker than the dif- 
ference in false-alarm probabilities, resulting in a partly 
more powerful detection statistic. To test this explana- 
tion, we have repeated the Monte-Carlo simulation with 
empty ranges (i.e. a fixed transient window), and with 
ranges where r^in > At", such that all waveforms over- 
lap in the transient start-time range, independently of 
window-function type. In both cases the "anomaly" dis- 
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FIG. 4: Detection efficiency for different transient-window 
types. Tlie piot sliows detection power of tlie Bayes factor B"^ 
and maximum-likeliliood statistic J-"max for assumed window- 
type vj in the case of injected signals of window type -cus. 
Transient parameters were drawn from range II of Eq. (671 
witfi a fixed SNR of p = 3. Tlie upper plot is for rectangu- 
lar injected transients (zu^ = r), while the lower plot is for 
exponential injected transients {wa = e). 



appears. Ranges I and II both have the feature that 
Tmin < At", such that the waveform-overlap will be sub- 
stantially lower for the rectangular-window compared to 
the exponential-window waveforms, resulting in a large 
difference in independent trials in the noise case. These 
observations are consistent with the explanation that the 
"anomaly" is caused by the difference in independent tri- 
als. Interestingly, however, this effect is not observed 
for the marginalized B ayes-factors, which always seem to 
correctly take into account the effective size of the pa- 
rameter space. 



C. Bayes factor self-consistency condition 

One can derive an interesting self-consistency condi- 
tion from the general definition (29) of the Bayes factor. 



which provides a useful Monte-Carlo test of our imple- 
mentation: the probability of obtaining a Bayes factor 
Bsg{x) e [Bq, Bq + dB] under any hypothesis Hi is 
given by the probability of obtaining a measurement x 



in the infinitesimal volume slice AVq = {x : Bsg{x) € 
[Bq, Bq + dB]} under that hypothesis, i.e. 



PiBo\H,) dB 



Pix\H^) d^x. 



(68) 



AVo 



Changing local coordinates from a; to y = {B{x),y±], 
where y± denotes n — 1 coordinates on the constant- Bq 
hypersurface So = {a; : -830(2;) = Bq], this can be 
written as 



P{Bq\H,) = f P{x\H{) dS, 



(69) 



where we defined the surface element dS = J d"^~^y± and 
we assume that the Jacobian J is non-singular, i.e. J = 
\dy/dx\ ^ everywhere inside AVq. Using the definition 
(29) of the Bayes factor, we can substitute P (xlHs) — 
j3sg{x) P (xIHg), and obtain 

P{BQ\ns)^ [ Bsg{x) P {x\nG) dS 

JSo 

= Bq [ p{x\nG) ds 

JSo 

= SoP(Bo|Hg), (70) 

where we used the fact that i3sG(3;)lsQ = Bq by defini- 
tion of Sq. We therefore find a general self-consistency 
relation for any Bayes factor, namely 



^SG(a;) 



P{x\Hs) PiBsGlUs) 



Pix\Hc 



P{Bsg\Hg) 



(71) 



In the framework of Monte-Carlo simulations [40], this 
implies that if we draw random data x according to 
the priors assumed in the Bayes factor, the ratio of 
probability-densities of obtaining i?sG = Bq in the signal- 
and the noise-case is identical to Bq. If the assumptions 
are satisfied, the signal-distribution of Bsg is therefore 
not independent of the noise-distribution, but is uniquely 
determined by it. If we know P {Bsg\T~I-g), then we also 
known the signal distribution, and vice versa. 

The resulting self-consistency relation for the odds ra- 
tio (pSl) is 



OSG 



PjOsGlHs) P{Hs\I) 

P{OsG\nG) pCHgII) ' 



(72) 



where the prior odds ratio determines the probability of 
drawing a sample x from the signal- or noise-population, 
respectively. Therefore the odds ratio predicts the ratio 
of event densities at any value OsC: rather than the ratio 
of normalized probability densities. 

Note that in order for Eq. ( 71 ) to hold for the transient- 
CW Bayes factor Bjr defined in Eq. (55), one must not 



draw signal amplitude parameters according to the 
physical priors (47), but according to the (unphysical) 
J^-statistic priors ( 54 ) that went into the construction of 
Bjr. The self-consistency relation (71) can equivalently 
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FIG. 5: Results of Monte-Ca rlo simulation of the Bayes-factor consistency relation (741, using A'^ = 10'' draws, transient 
range II, and amplitude priors (541 with cutoff pmax = 7 {left panel), and pmax = 14 {right panel), respectively. The upper plots 
show the distributions of logBjr in the noise-case (dashed), and in the signal case (solid), and the lower plots show log(ps/PG), 
which should coincide with logBjr (dot-dashes line) according to (|74[). We show the binned histogram values (thin lines) as 



well as kernel-smoothed fits (thick lines), and we have restricted the consistency test to a region of good overlap (indicated by 
vertical dashed lines) between both distributions in order to avoid numerical problems. 



be expressed as 
logB^ = logP(logB^|Hs)-logP(logB^|HG) , (73) 

which is more directly suitable for numerically testing 
this relation in a Monte-Carlo simulation. Defining the 
shortcut Pi = P (logSjrI'Hi), this can also be written as 



log 



PS 
Pg 



logS^. 



(74) 



We have performed a Monte-Carlo simulation generating 
values of Bjr in the noise- and signal-cases, with ampli- 
tude parameters drawn according to the (unphysi- 
cal) J^-statistic priors (54). Figure [s] shows the resulting 
distributions of pg and pq, and the plots of log(ps/pG) 
versus log Bjr, which should fall on a straight line of 
unit slope according to the self-consistency relation (74). 



These results show that the self-consistency relation is in- 
creasingly well satisfied with increasing prior cutoff Pmaxj 
in particular we find good agreement for cutoff values 
above Pmax ^ 10, as illustrated in Fig.js] This can be un- 
derstood as follows: for smaller values of Pmax, the noise 
population P {BjtIT-Lq) is biased towards larger values of 
Bjr, because the approximation in Eq. (49 ) is increasingly 



violated. In the noise-only case, the likelihood ratio C 
will peak somewhere around = and fall off accord- 



ing to a Gaussian (31 ) with characteristic width of order 
cr(p) ~ 0{1), modulo geometric factors or order unity. 
Therefore the value of £ will not be negligible at the 
cutoff boundary Pmax ^ (l)i and the extension to in- 
finity will overestimate the integral. Therefore Pmax S> 1 
is necessary for the integral to be well approximated by 
Eq. (pgl). 



VI. PARAMETER ESTIMATION 

Parameter estimation simply consists of computing the 
posterior probability for the signal parameters 9, given 
the observed data x, namely 

p{9\x,ns) - cP{x\ns,0) p{0\ns) , (75) 

where c = l/P{x\'Hs) is a normalization constant inde- 
pendent of 6. Often one is not interested in a simulta- 
neous estimate of the full set of parameters 0, but only 
in a subset 9i C 9 without regard for the "nuisance" 
parameters 02, where 9 = {9i,92}- From the general 



14 




5 10 15 20 25 30 




5 10 15 20 25 30 





0.3 




0.25 


m 


0.2 






8 


0.15 






aT 


0.1 




0.05 








Oh 



- 




III- 
Po = 5 


— / 




1 \ 



2 4 6 8 10 12 14 




FIG. 6; Example posteriors (791 on t" (left column) and r (right column), for an injected rectangular transient-CW signal within 
range III, with randomly drawn amplitude parameters at fixed optimal SNR po = 5 (upper row) and Po = 8 (lower row). The 
solid vertical line indicates the injected parameter value, and the dashed vertical line ('ML') indicates the maximum-likelihood 
estimate. 



expression P {6i\. . .) — J P {9i,62\ ■ ■ ■) (102 and Bayes' 
theorem (|75[), we obtain the marginahzed posterior 



p{ei\x,ns) cc j P{x\n^,ei,e2) p(0i,02|Hs) de2. 

(76) 

For our present model the hkelihood function ( 26 ) can be 
expressed as P {x\9) = ne'^^^^^l"^ C(x] 9) in terms of the 
likelihood ratio C of Eq. (31). Using J^-statistic priors 



(54), we can perform the yl^-integration explicitly and 



obtain 



P (T, A|a;, Hs) « e^(-^^^-^) P [T , A|Hs) 



(77) 



which is a useful starting point for further marginaliza- 
tion. If we consider a targeted search in Doppler param- 
eter, i.e. P (AjHs) = ^(A — As), with an assumed window 
function type w, we can write the posterior probability 
for r} as 



P {t", t\x, Hs, a, vj) cx e^(=^'^''^' P {t\ t\Hs 



(78) 



and the respective marginal posteriors on the transient 
parameters are simply 



/'(t°|a;,-Hs,A,n7) cx J e^^'''^''^^ dr , 
P(T|a;,Hs,A,w) (X J e^'^'^'^'^Ut" , 



(79) 



where we assumed uniform priors (46) for {t'^,T}. The 



generalization to marginalization over the window type 
m is straightforward, and yields a weighted sum of these 
posteriors with relative prior probabilities of the different 
window types, e.g. 

P {t°\x, Hs, A) cx ^ P {t''\x, Hs, A, vu) P (n7|Hs) . 

m 

(80) 

Similarly, parameter-estimation on the window-type vj 
itself can be expressed as 



P {vj\x, -Hs, A) cx P (ro|-Hs) J e-^(^'^''^) dt° dr 
cx P(n7|'Hs) Bjr{x;X,w), 



(81) 



so the window-type specific Bayes factor ( |55[ ) is propor- 
tional to the relative likelihood of different window-types 

TU. 

In the frequentist framework one often uses maximum- 
likelihood estimators for parameter estimation, i.e. 
{^ML' ■^ml} such that 

•^(a;;^ML>^ML) = max J"(a;;i°,T) , (82) 
for fixed Doppler-point A and window-type w. 
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FIG. 7: Parameter-estimation accuracy: rms-errors on t° {left panel) and t {right panel) as a function of SNR for maximum- 
likelihood (ML) and maximum-posterior (MP) estimates. Each value of SNR corresponds to TV = 10^ random parameter draws, 
and the error-bars represent la jackknife error-estimates. The dotted horizontal line ('rand') corresponds to a uniform random 
guess (851, and the dot-dashed horizontal line corresponds to a maximally biased r = T^in "guess" (86 1. 



The following Monte-Carlo studies use rectangular 
transient windows within a wider range than ( 67 1 , 
namely 

range III: re [0.5, 14.5] days , t° € [0, 30] days . 

(83) 

Using the transient-parameter range III, Fig. [6] shows 
one example of parameter posteriors ( 79 1 and maximum- 
likelihood (ML) estimators (82) on t'^ and r for injected 
signals with SNR — 5 and po — 8, respectively. We 
see in Fig. [g] that the timescale of variations in P {t^\x^ is 
shorter than in P{t\x). This is due to the combined ef- 
fect of the wider plotted range in t*^ and the twice smaller 
characteristic correlation timescale in t", as derived in 



Eq. (63) 



In order to study the quality of transient parameter- 
estimation as a function of SNR, we performed Monte- 
Carlo simulations comparing maximum-posterior estima- 
tors (MP), defined as 



P{t°Mp\x,ns) =m^axF(t°|a;,Hs) 
P(rMp|a;,'Hs) = maxP(T|a;,'Hs) 



(84) 



and maximum-likelihood estimators (ML) of Eq. ( 82 1 , 
using uniform priors on {t^ , r} within range III, and 
physical priors (47) on amplitude parameters. For dif- 
ferent fixed values of SNR we perform TV = 10^ simu- 
lated parameter estimates, and compute the rms errors 
from i° — and t — Tg, for ML- and MP- estimators, 
respectively. We also compare the results to the error of 
a pure random guess within the range [ti , ti + At] , where 
At = t2 ^ ti. For uniform priors this is 



(rms [t-ts])' = 



1 

Ai2 



dts / dt{t-ts)^' 



At" 



(85) 

For the transient range III this yields random-guess er- 
rors of rms \iP — tf\ ~ 12.25 days, and rms [r — Ts] w 



5.72 days, which are shown in Fig. [T] On the other hand, 
for a maximally biased "guess" of r = r^-ain i one finds 



rms [rniin — Ts] = At/ a/S ~ 8.08 days . 



(86) 



The results of the parameter-estimation Monte-Carlo 
simulation are shown in Fig.[7] For low SNRs of p < 7, we 
see that the MP-estimators perform better than the ML 
estimators, while for higher SNR the estimation quality 
of both estimators converges. We note that for p -> the 
parameter estimation on t° converges to a random guess 
as expected, but in the case of tml we notice a substantial 
deviation, and to a lesser degree, also for tmp- These esti- 
mates fall closer to a maximally biased r = T^-m "guess" , 
which indicates an increasing bias in the t estimators for 
low SNR, strongly favoring values close to Tmin- This 
surprising effect will be studied in some more detail in 
the following section. 



A. Estimation bias on timescale r in pure noise 

Figure [8] shows normalized histograms of the parame- 
ter estimates on r in pure Gaussian noise, i.e. for p = 0. 
These results confirm the estimation bias towards Tmin 
previously seen in Fig. [7] We have been able to trace this 
effect to a surprising fundamental feature of the Gaus- 
sian random walk underlying matched filtering. If we 
discretize the integration time as Tj = j AT in steps of 
AT, then linearity of the scalar product (32 1 implies 



<(?;■) = <(7;-i) + Ax' 



(87) 



where Ax'^ ^ is independent of x'^{Tj_i) and follows a 
Gaussian distribution with zero mean (in the noise case 
Hg)- We see from this expression that the amplitudes x'^ 
can be interpreted as Gaussian random walks over finite 
steps AT in the integration time. The four Gaussian 
random walks {x'^} are combined in the quadratic form 
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FIG. 8: Parameter-estimation bias in tml {left panel) and 
tmp {right panel) in pure noise data (SNR=0). The plots 
show normahzed histograms of A'' = 10* parameter estimates 
of T in Gaussian random noise data. 
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FIG. 9: Distribution of random-walk index n„ 
maximum of occurs (A*' — 10* trials) . 



100 



at which the 



(361 to yield T{Tj), where they are normalized by A^^"^ 
such that J- follows a ^^-distribution with four degrees 
of freedom for any Tj , independently of the random- walk 
step j. 

We can therefore consider a simpler toy model, namely 
a 1-dimensional normalized random walk, defined as 



1 



2. Ax,. 



where the Axj are n independent Gaussian random vari- 
ables with zero mean and unit variance, i.e. Aa;, ^ 



Gauss(0, 1). The random walk (88) is normalized in such 
a way that it follows exactly the same distribution at ev- 
ery step n, i.e. s„ ~ Gauss(0, 1). We denote as s"^^^^ the 
maximum of s'^ over a walk-step window [rii, 712], i.e. 



max , 

nS [rii ,n2] 



(89) 



and plot the distribution of rimax & [ni , n2] in repeated 
trials of such normalized random walks. For example, if 
we consider ni — 50 and n2 — 100, and TV = 10* trials, 
we obtain the distribution of maxima shown in Fig. [9j 
which illustrates a qualitative bias towards rimin, and to 
a lesser extent n,„ax, similar to what was seen in Fig. [8] 
for the physical parameter estimation of r. This seems to 
be a manifestation of a well-known property of random 
walks, namely Levy's arcsine law, which is discussed, for 
example, in Sec. 4. 2 of [4T] . 



VII. CONCLUSIONS 

We have discussed the case for a transient-CW search, 
from the point of view of its astrophysical motivation, 
and in order to bridge the gap between short burst-like 
signals and traditional infinite-duration CW signals. 



We have introduced a simple transient signal model 
based on the classical CW signal, modulated by a window 
function of finite support. The corresponding Bayes fac- 
tor has been derived and implemented, and we have per- 
formed Monte-Carlo simulations to compare its efficiency 
to the orthodox maximum-likelihood detection method. 
These results show that the Bayes factor is both more 
sensitive and more robust than a maximum-likelihood 
statistic, and it yields better parameter estimates, at sim- 
ilar computing cost. 

The Monte-Carlo studies presented in this work have 
been limited to pure Gaussian noise, and it would be im- 
portant to test this method in practice on actual detector 
data. In particular one needs to address the question of 
how to assign meaningful false-alarm probabilities to de- 
tection candidates in real detector data (see also [5 for 
an example of these difficulties). 

The search-method discussed here is currently re- 
stricted to non-repeating transient CWs, and more work 
would be required to generalize this approach to allow 
for repeating transient-CW signals. 

In addition to the fully coherent search method, we 
have derived the necessary formalism for a semi-coherent 
transient search, which could be used to perform an all- 
sky, all-frequency wide parameter-space transient search, 
for example running on Einstein@IIome. More work is 
required to fully develop and implement this approach. 
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Appendix A: Transient search Implementation 

Our numerical implementation of the Bayes factor 
Bjr{x; A) of Eq. (55) consists of two steps: 



1. calculate a discretized J^-statistic map J^{x;t'^,T) 
over the search ranges in t'^ and r, 

2. compute Bjr{x) by discrctizing the marginalization 
integral in Eq. (55) as a sum. 



The following two sections provide more details about 
these two steps, respectively. 



1. Atoms-based J^-statistic computation 

We discretize the 2-dimensional J^-statistic map over 
the search ranges f G [^min'^min + ^^"^l ^^"^ ^ 
[Tmin , Tmin + At] in stcps dt'^ and dr, respectively. 
Namely, we compute the J-'-matrix 



over the Nfo x Nt rectangular grid 



''m ''mm "i- ■ 

Tfi — Tinin ~t~ ^ dv , 



(Al) 



(A2) 



where Nto = At°/dt° and = Ar/dr. 

The default for these step-sizes used here are dt'^ — 
dr — 1800 s. If the shortest signals considered are 
Tmin ~ 0.5 days long, then according to Eq. (62) the 
worst-case mismatch is m ^ dr/r^in ~ 0.04, i.e. a 4% 
loss of squared SNR. 

In the current implementation of the transient- 
CW search we use the underlying discretization of 
the J^-statistic computation in ComputeFStatistic_v2 
(CFSv2), as described in more detail in [20]. Namely 
the J"-statistic (36) is computed from M'^^ and x'^ of 
Eqs. (32), (33), which are approximated as sums 



iVsFT 
iVsFT 

1=1 



(A3) 



in terms of the J^-statistic "atoms" 

f ti+TsFT 



M 



X 

■E 

X 



x''{t)hf^{t)dt 



(A4) 



h'^it) h^{t) dt, 



where Tsft is the length of the short Fourier trans- 
forms (SFTs) that are used as input data, typically 
Tsft = 1800 s. In the above expressions we implicitly 



assumed that the transient-window function g{t) varies 
slowly and can be approximated as constant over the 
timescale Tsft- For any chosen Doppler position A, the 
code first computes the TVsft atoms {a;^^^, ^ 
over the whole observation time of interest, which are also 
the primary input to this implementation of the standard 
CW J^-statistic. The J^-statistic value Fmn for any par- 
ticular transient parameters {t^, t„} is then computed 
from the corresponding partial sums in Eq. ( A3 ) . This 



approach allows for an efficient computation of the J-m n 
map in the case of a rectangular window function of 
Eq. (17): going from t„ to t„+i can be achieved by a 



single extra addition, namely 



(A5) 



(and similarly for Ai^u), where ii is the atom-index cor- 
responding to the time-step + t„ (assuming for sim- 
plicity that dr — Tsft)- 

In the case of the exponential transient window gc of 
Eq. (18), the whole sum in Eq. (A3) needs to be re- 



computed for every matrix-element m, n, as the window- 
function provides different weights at every point. 



2. Transient marginalization integrals 



Given the transient matrix J- mm we can now turn the 



integrals in Eqs. ( 55 ) and ( 79 ) into simple sums. A minor 



subtlety arises because of the potential numerical prob- 
lems of expressions like , which overfiow in double pre- 
cision for values F > 709. Such values are easily possi- 
ble for noisy non-Gaussian data with line-artifacts or for 
strong injected signals. It is therefore numerically safer 
to rewrite these sums using the discretized maximum- 
likelihood statistic value ( 53 ) , namely 



(A6) 



and write Bjr of Eq. ( 55 ) as a discretized sum in the form 
\ogBjr{x; A) « J'm 



Co 



E E ' 



, (A7) 



where we defined cq = log(70/(pj^ax -^t" ^t)), and 

A^ }-rm — J'mn max — . (-A^) 

For large values of J^raax, some terms e^"^™" can now 
underflow to zero, but this poses no numerical problems 
because these contributions were negligible anyway. Sim- 



ilarly, we can compute the parameter posteriors ( 79 ) as 



cx 



n=l 



P{Tn\x,'Hs) OC E 



(A9) 



(AlO) 



m—l 



where J^max only affects the normalization and has been 
dropped. 
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3. Computing cost 

In order to be able to plan which types of searches can 
be performed with a reasonable investment of computing 
cost, it is useful to have a rough computing-cost model 
that allows one to predict the expected run time of any 
search. 

Let us start with the underlying J^-statistic imple- 
mentation CFSv2 which is used in our current coherent 



depend on the start-time index m. Therefore we obtain 
the total J^-map cost as 



J-'map 



At° At (r^in + At/2) 



dt° dT 



Ts 



FT 



(A13) 



transient-search implementation described in Sec. A 1 



which is quadratic in At (and where we assumed 

Nr > 1). 

In the case of a rectangular transient-window, various 
sums are closely related and we can use the optimization 



Note that more efficient J^-statistic algorithms do exist, 
based on resampling and FFT techniques |321 SSI ■ How- 
ever, the "Williams-Schutz" method currently im- 
plemented in CFSv2 is well suited to our purpose, as it is 
already based on computing the atoms Al^i/^i} over 
SFTs. Nevertheless it will be interesting to study possi- 
bly more efficient transient-CW implementations based 
on resampling and FFT techniques. 

We can break up the total computing time c^' for 
Bjr{x] A) at one Doppler point A as follows: 



of Eq. ( A5 1 to reduce the computing cost. Namely, every 
step n — > n -I- 1 in the timescale r adds just the cost 
of one extra time-step dr. This cost is Cr dr, where Cr 
is the machine-dependent time to do all sums for a unit 
timescale in the case of a rectangular window. For every 
start-time index m, we need to compute the sums up to 
''"min first, costing CrTmin. Summarizing, we can express 
the accumulated computing cost c5„ = J2n n P^r line 



of the J^m n matrix as = Cr Tmi, 
The total J^-map cost is therefore 



Cr NrdT — Cr T„ 



1. the time Catoms to compute the atoms (A4|, 



c = c 

J-'map 1' 



At" (i 

IF ^ 



At) 



TSFT 



(AM) 



Bayes-factor marginalization: the marginalization 



2. the time c^^^^g^p to compute the J'-map (Al) over 
{i",T} for given window-type nj, and 

3. the time Cmarg to marginalize over the J^-map to ^nd so we can directly write the marginalization cost 



( A7 ) is a simple sum over the exponentiated J^-map ma- 



obtain the Bayes factor Bjr{x; A), using Eq. (A7) 



Atoms cost: the time Catoms to compute for one 
Doppler position A the A^'sft — Tdata/?SFT atoms using 
CFSv2 is simply 



Cl 



At° At 
W ~d^ 



(A15) 



Catoms — ^0 



^ data 



FT 



(All) 



where Cq is a machine-dependent timing constant. 

jF-map cost: the time c^j^^^p to compute the Nto x Nr 
matrix of values {Fmn}^ where iVjo = Ai"/dt°, and Nj- = 
At/c^t. This time can be expressed as 



where Ci is the machine-dependent cost of exponentia- 
tion and summation of real numbers. Note that expo- 
nentiation is a very costly operation, and so we use a 
lookup-table approximation to reduce this cost. 

The total computing cost c^'' for the transient Bayes 
factor at one Doppler point A is now expressible as 



^atoms i '•JF'map ' 



(A16) 



f 

^J-'map 



EE' 



(A12) 



and for a search over N\ Doppler points, this would sim- 
ply extend as 



= iVAC 



(1) 



where c^ „ is the time required to compute Fm n for par- 
ticular values m, n, which depends crucially on the type 
tu of window-function. 

Let us first consider the exponential transient-window 



(A17) 



(18 1, which corresponds to the generic case where no 
special optimizations can be used. Note that at the 
smallest timescale index, n = 1, we already had to 
sum all atoms corresponding to a timescale Tmin. How- 
ever, for every matrix element we need to re-compute 
this sum, as the window-weights will be different every 
time. This introduces a constant computing-time offset 
c„^g — Cc Tinin, where Cc is the machine-dependent time to 
compute all atom-additions and window-weights within 
one unit-time. From this we can express the cost for com- 



where the extra cost of summing these Bayes-factors will 
be negligible. 

We have verified that these timing models are a good 
description of the actual performance of the code [48] 
by varying the parameters iVspT, At" and At, and by 
fitting the measured times to the model. This yields the 
following timing constants on an Intel Core2 Duo CPU 
with 2.60 GHz (Lenovo T61p): 



Co = 1.4 X 10"" sec, Cl = 2.8 x 10"" sec . 



Cr = 4.2 X 10"* sec, Co = 1.3 x 10"' sec 



-7 



(A18) 



puting J"„„ as c° 



Co (Tmin + TL dT) ^ which docs not 



If we target a known pulsar with a transient-CW search 
using data from 2 detectors spanning one year (A^sft = 
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35,000), with a timescale range of r € [0.5, 14.5] days, 
using step-sizes dt° = dr = Tsft = 1800 sec, we obtain 
an estimated computing costs of 



Catoms ~ 0.05 sec, Cmarg ~ 0.3 SCC 

c>map ~ 0-5 sec, c5r„,ap « 540 sec . 



(A19) 



We see that such a targeted search would be easy to 
perform for all interesting pulsars, even using the much 
slower transient exponential window. Comparing this 
transient-CW search to a coherent CW search over one 
year, we notice that the cost per template A is about 18 
times higher than the CW search for a rectangular win- 
dow, and about 10^ times higher for the exponential win- 
dow. Wide parameter-space transient-CW searches will 
therefore be severely limited by computing resources, and 
a semi-coherent transient search method as discussed in 
Sec. |IV G| will be required. 



4. Synthesizing Monte-Carlo draws 

Following the method used in [34] . we have imple- 
mented an efficient Monte-Carlo simulation method, by 
avoiding the generation of the primary data-input of the 



search code (time-series or short Fourier transforms), and 
instead synthesizing higher-level secondary data-input to 
the transient-search functions directly. In our case, we 
synthesize the iVsFT atoms {x^^^i, Mfi^^i}, which arc the 
intermediate input-data to the transient-search function. 
This approach is very economical in computing resources 
and allows us to generate large numbers of Monte-Carlo 
draws in a very short time on a single machine. We 
draw signal parameters {^s, As, 7^} according to the pri- 
ors, and from this we can compute the deterministic 
signal-atoms s^^i of E g. (| 38[) and the antenna-pattern 
atoms M 




The noise-atoms n^^i are 



t^^,j of Eq. ( |A4 
Gaussian random variables with zero mean and covari 
ance matrix M 



These can be generated from un- 
correlated Gaussian variates, e.g. by using a Cholesky 
decomposition on A4p,i,^i. The data-atoms are then sim- 
ply i ~ n^^i + Sp^i, according to Eq. (37). Note that 



we have assumed the transient-window function to be 
constant on the atoms timescale Tsft, and therefore we 
can synthesize standard non-transient CW atoms. The 
transient-window function is applied when computing the 
J^-statistic map, cf. Sec. |A 1[ All Monte-Carlo simula- 
tions used the "Mersenne Twister" random-number gen- 
erator gsl_rng_tiitl9937 from GSLfiS]. 
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